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Studying charge current through weak links, the 
Josephson effect, is one of the most important-gart of su- 
perconductivity, both theory and experimentuD. Besides 
general interest, this problem is important for engineer- 
ing the numerous devices based on the Josephson effect. 
The Josephson effect reveals itself in tunneling junctions 
as well as more complex mesoscopic structures built of 
superconducting and normal layers. A major part of 
the theoretical results in this field has been obtained |USj- 
ing the method of quasiclassical Green's functionscloiJ'ElLl. 
The advantage of this general method is that disorder 
and inelastic processes can be conveniently incorporated 
into the theory. Iii-a ballistic case, one can apply a more 
simple techniqueto based on the Bogoliubov - de Gennes 
(BdG) equationty. A strong side of the BdG-approach is 
that it is valid for description of the interface reflection 
and transmission, where the potential varies on a micro- 
scopic length and the quasiclassical .-theory in its origi- 
nal form fails. As shown by ZaitsevM, an isolated par- 
tially transparent interface can be taken into account by a 
proper boundary condition for the quasiclassical Green's 
functions. Recently, there has been considerable techni- 
cal progress where the boundary conditioH-is formulated 
using the Schopohl-Maki paramcterizationlla of the quasi- 
classical Gi-pen's functionEJ, or in terms of effectiiia wave 
functionstllEj. However, it has been arguedllsllZI that 
these boundary conditions may give a wrong result in 
the case when a coherent scatteiing by several interfaces 
takes place. It is shown in Ref.tj that the quasiclassical 
density of states of a SS' sandwich disagrees with the 
"exact" one, found from the Gor'kov aquation. The dis- 
agreement has been attributed in Ref.EZi to the presence 
of closed trajectories formed by sequential reflections by 
the interface and the outer boundaries of the sandwich. 
In particular, the correction to the quasiclassical Green's 
function due to the loop-like trajectories violates the nor- 
malization condition, which is an essential element of the 
quasiclassical technique. In this paper, we extend these 
results to the case of an open geometry and analyse ap- 
plicability of the quasiclassical theory to the Josephson 
effect in a multi-layer mesoscopic structure. 

The quasiclassical theory is a simplified version of 
the "exact" theory of superconductivity based on the 



Gor'kov Green's function formalism. The main assump- 
tion made in the course of its derivation is that the poten- 
tials vary slowly on the Fermi wave length Xp = '2t^/pf, 
Pf being the Fermi momentum, that is the parame- 
ter Xp/^Q is small, where is the coherence length 
(Co ~ vp/^, vp is the Fermi velocity). The question we 
address in this paper is whether there are corrections to 
the theory which are not controlled by the quasiclassical 
parameter Xf/£,o ^ 1. 

To judge if the quasiclassical approach gives valid re- 
sults, we compare its predictions with the solution to the 
BdG equation. In the clean case of the mean-field theory, 
the BdG approach is fully equivalent to the Green's func- 
tion method, which is the starting point to the derivation 
of the quasiclassical approximation. For this reason we 
consider the BdG solutions as "exact" for the purpose of 
the comparison. More specifically, we consider a multi- 
layer SS'S" . . . structure shown in Fig.|| and calculate the 
angle and energy resolved partial current j{0,e), where 
6 is the angle of incidence and e is the energy of the ex- 
citation propagating through the multi-layer structure. 
Comparing the results of the two approaches, we make 
our conclusions on the validity of the quasiclassical ap- 
proximation. 

As discussed in detail in Ref.0lli, one cannot make 
any conclusions comparing directly jBdoiOi evaluated 
from the BdG-equation with its quasiclassical counter- 
part jqc{0,e). The point is that in the BdG approach, 
the incident particle is taken as a plane wave with pre- 
cisely defined wave vector whereas in the quasiclassical 
approach one deals with classical trajectories, where the 
momentum in the direction perpendicular to the velocity 
has quantum uncertainty. The infinitely extended plane 
wave suffers multiple refiections on the interfaces, and 
the refiected/transmitted waves inevitably interfere be- 
cause of the infinite extension. The interference leads to 
an intricate picture of Fabry-Perot like resonances and 
a fine structure in the angle dependence on the scale 
69 ^ Xp/a where a is the layer width. To illustrate 
this point, we show in Fig. ^ angle-energy resolved jBdG 
current through a SISIS-system of superconductors (S) 
separated by two barriers (I) in a narrow region of angles 
(the BdG calculations are done by the method presented 
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FIG. 1: A planar, multi-layer structure consisting of A'' lay- 
ers separated by barriers and sandwiched by half-infinite elec- 
trodes. The complex order parameter denoted Ao, . . . , Ajv+i 
are considered as inputs. It is assumed that each layer is 
connected to an independent current source Io,Ii ■ ■ ■ so that 
one achieves any given distribution of the phase of the order 
parameters without violation of the current conservation . 
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FIG. 2: Angle resolved current in a SISIS structure as a func- 
tion of angle 6. The energy e — 1.2Ao. The order parameter 
in the leftmost superconductor is Ao. The order parameters in 
the next layers read Ai = e"^'''*Ao, A2 = iAq. The thickness 
of the internal layer is a = uf/Aq. The interface transparen- 
cies Ti = 0.1, r2 — 0.5. The dashed and solid lines show the 
BdG jBdG and the quasiclassical current jqc, respectively. 



below, see alsoEfl). The "exact" current shows rapid and 
strong fluctuations in the region where the quasiclassical 
current is almost a constant. However, on a large scale of 
angles, judC averaged in a small angle window (ca 
grain current) is a smooth function. As discussed intuit 
the coarse-grain averaging is equivalent to building sta- 
tionary wave packets, peaked on classical trajectories, on 
which the quasiclassical theory is formulated. It is on 
this low resolution level where the quasiclassical and ex- 
act theory are expected to agree with each other. For 
these reasons, we use only coarse-grain BdG-current for 
comparison with the quasiclassical theory. For definite- 
ness, we calculate the current at the leftmost interface of 
the multi-layer structure. 

In principle, the coarse-grain averaging can be per- 
formed analytiaajly applying the path length expansion 
method of Ref.tZl. (In case of a double-layer, the aver- 
aging can be done directly with the result expressed via 



the elliptic integralst^JHU.) Nevertheless, we perform the 
averaging by a numerical integration to avoid lengthy al- 
gebra. 

The paper is organized as follows. In Sect.| we dis- 
cuss solutions to the Bogolubov - de Gennes equation 
In Sect 



lA 



we build plane waves and then, in Sect.IB, 



we introduce the scattering S-matrix. In Sect. I C, we ex- 
press the BdG current via the elements of the S-matrix. 
The method which we use to evaluate the current in the 
quasiclassical technique is presented in Sect.||. Numeri- 
cal results are shown in Sect. |l|, and their interpretation 
is presented in Sect.IV. In Sect.^ we discuss validity of 
our results in more realistic models. Technical details of 
the derivation are collected in Appendices. 



I. BOGOLUBOV-DE GENNES EQUATION 

In this section, we consider the theory of multi- 
layer structure in the framework of the Bogoliubov - de 
Gennes equation. Stationary two-component wave func- 
tion ifiir) = (") of an excitation with the energy E sat- 
isfies the BdG-equatioiS, Hip — Eip, 



H = 



f ^{p-^A) + V 
\ A* 



A 



(1.1) 



1^, pf being the Fermi momentum. 



where ^{p) = §^ 

V{r) is the potential energy, A(r) is the complex order 
parameter, and A{r) is the magnetic vector potential. 
The charge current density J can be found as 



J 



where J 



-c^ is the current operator 



Tz-A 

c 



(1.2) 



(1.3) 



Tz being the Pauli matrix. 

The non-diagonal elements of current operator J(r), 
{'ipn\J{r)\'ipn') are evaluated as 



Jlpn tpn' + 



(1.4) 



In superconductors, the charge current created by an el- 
ementary excitation in a state ipn is not a conserving 
quantity, i.e. divJ„„ 7^ 0. The charge conservation is 
restored after the summation over all the BdG excita- 
tions, provided the pair potential A is self-consistent. 

To take advantage of the unitarity property, one con- 
siders the conserving quasiparticle current, j'^^ir), cal- 

that 



culated with the help of the operator, j'^p = 



— (tzP - -A) 

TO V C / 



(1.5) 
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and 

The continuity equation 

div = . 



(1.6) 
(1.7) 



follows from the BdG equation. 

We solve the BdG equation for the case of a planar 
structure shown in Fig.|l|. It is composed of N layers 
surrounded by two half-infinite homogeneous supercon- 
ductors. The complex order parameters in each layer is 
a constant which is taken as an independent input. 

Choosing the x axis perpendicular to the layer plane, 
the partially transparent interfaces are modelled by the 
(5-function barrier 



V 



A 



5{x) 



(1.8) 



where A is the strength of the potential, m is the mass. 
Below, we characterize the interface by the parame- 
ters R = \\/{pF + «A)|2 and T {R + T ^ 1) which 
have the meaning of the reflection coefficient (R), and 
transparency(r) for the normal incidence. 



A. Plane wave solutions 

Due to the in-plane translational invariance, the solu- 
tions can be taken in the following form 



(1.9) 



where p|| is the in-plane momentum, and ^ obeys the 
one-dimensional BdG equation. 

First, we consider the plane wave solution to the BdG 
equation in each of the layers where A = const and V — 
0. The hmction ip{^) satisfies the equation 



Tpix) — Eip{x) 



(1.10) 



where = {px+Pf\^PF^ f^m, and = —id/dx. 



Eq.(l.lO) has 4 linearly independent plane wave solu- 
tions: 

^,^{x)=e"'P^'=^, , u = ±, a = ± (1.11) 
where the momentum pi, is found from (see Fig. ^ 

P. = ^pI - pI ± 2m^ , 3? p± > (1.12) 

with 

e= Vi?2_|A|2 . (1.13) 

We choose the branch of the square root so that ^ > at 
E > |A| and 9^ > when E < |A| 




FIG. 3: 4 linearly independent solutions to BdG equation. 
Index V describes the type of the quasiparticle v — + is the 
electron-like and u = — the hole-like excitation; a defines the 
sign of the a;— component of the momentum. As shown by 
arrows, the excitation propagates to the right if i/ ■ a = +1, 
and to the left otherwise. 



The two-component amplitudes ip± found from 



Eq. ( 1.10 ) with substituted for may be chosen in 
the following form 



c 



where 



2^ 



E + (, 



(1.14) 



(1.15) 



(c = \/l — ab). These expressions are applicable for any 
energy E including the gap region E < |A|. Outside the 
gap, b = a*. 

The amplitudes are normalized to the unit flux, i.e. 



W±Tzi^± = ±1 



^> |A| 



(1.16) 



Note the orthogonality relation (outside the gap), 
'0±'rz'0=F ~ 0, in agreement with the current conservation 
in Eq. ( |1.7| ). For any E, the determinant of the matrix 
[■(/'+, V'-J the columns of which are V'+ and V'-, equals to 
unity (in other words, ip'^iTyip- — 1). 

For future needs, we define conjugated amplitudes 



(1.17) 



which posses useful properties of orthonormality and 
completeness: 



1 



(1.18) 



where 1 is the unit 2x2 matrix. Outside the gap, where 
4>l = tpt Tz , the orthonormality expresses the quasipar- 
ticle current conservation. 
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The physical meaning of the quantum numbers v and 
cr is clear : in a ccordance with the sign of the probability 
flux Eq. (1.16), the excitations with = 1 are electron- 
like whereas the states v — —X are hole-like. The pa- 
rameter a shows the direction of the momentum. The 
excitations for which the product av equals to +1 (- 
1) propagate in the positive (negative) direction of the 
x— axis. 



B. The scattering matrix 



J(x) being the current operator Eq. ( |l.3| ). We restrict 
ourself to the simplest case where the distribution func- 
tion depends only on energy, i.e. ni{E) = n{E). Then, 
the current can be written as 

oc 

J{x)^ j dp\\ j dE{2n{E)-l)J{E,p\\;x) (1.24) 



where the partial current density, J(£',p||; x) at the point 

X, is 



To introduce the scattering matrix, we first define the 
in-coming and out-going free states, the amplitudes of 
which are related by the S-matrix. We number the in- 
coming plane wave states in the following way: 



|3) 



(0) 



As before, the first of the lower indices of ^''s specifies the 
electron- hole degree of freedom ( "-)-" for electron, and "- 
" for hole) and the second one shows the direction of the 
momentum a — ±; the upper index L or R specifies the 
initial location of the excitation on the left or right side 
of the structure. The out-going states are 



1 3) out 



- + 
{R) 



The basis wave function read 



(1.20) 



(1.21) 



with ■01^ from Eq. (1.14), these function arc normalised 
to 5{p\\-p\^)5{E-E'). 

In the L- or R-regions of free motion, the solution to the 
BdG equation, |i), corresponding to the incident quasi- 

particle in the state can be presented as 

4 



(0) 
out 



(1.22) 



These equations with i ~ 1, ... ,4 define the 4 x 4 S- 
matrix. The method which allows us to evaluate the 
elements of S-matrix is presented in Appendix |^ 

C. The BdG current 

Expressed via the distribution function of the excita- 
tions ni{E), ( where E is the energy and i — 1, . . . , 4 and 
E is the quantum number introduced in Eq. ( |1.22D ) , the 
charge current in the x— direction readsE3 

oo ^ 

J{x)= J dp\\ j dEY,{2n,{E)~l){i,E\J{x)\i,E) , 

(1.23) 



J{E,p\\-x) ^Y.^i,E\J{x)\i,E) (1.25) 

i=l 

To evaluate the current i n the left or right electrodes, 
we substitute \i) from Eq. (1.22), and take into account 



(1-19) the unitarity property, S*f,^Sfi = (5///. We get 



J(i?,P||;a;)^ jf (x) + 25R ^ S},Jj'> {x) (1.26) 



4 



where J 



(0)/ 



E (^) ~ Jkk{x): here summation is per- 

fc=i 

formed over the 4 plane wave states Fig. g on the left or 
right side of the structure. The meaning of is that it 
would give the (partial) current in the left or right region 
if the plane wave states with the given energy E were 
equally occupied. This is the contribution to the current 
which produces the bulk supercurrent 2eNsVs. In our 
case, Je = since the phase of the order parameter is 
assumed to be a constant within the outside regio ns. 

The second term in the right hand side of Eq. ( 1.26| ) 
is due to the interference of the incoming and outgoing 
waves. Considering for definiteness the left region, the 
initial states i = 1,2 interfere with the final / = 3,4 
states. For the energy E outside the gap, jj-^'' is other 
than zero only ifi = l,/ = 2ori = 2, / =:1 (i.e. for the 
interference with the Andreev reflected particle). The 
partial current density J{E,p\\) at the point adjacent to 
the first interface, i.e. at a; = 0~ reads 

J(i?, Pll) = 25R (52*iJi;^(0-) + Sl,jf^{0-)) . (1.27) 



Calculating the curren t mat rix ele me nts Eq. (A3) with 
the a mpli tudes in Eqs. ( 1.14 ), and ([A^), one derives from 
Eq. (|1.27D that 



J(ii;,P||)==-3?i(52iAo-^i2A5) , i?>|Ao|. 

(1.28) 

where Aq iSj-the order parameter in the left electrode 
(see Fig. |])lj. The scattering matrix is calculated by 
the tra nsfer matrix method as described in Section 

Eq. ( 1.28| ) gives the BdG current carried by plane wave 
states with definite value of p| | . This quantity strongly 



5 



fluctuates (see e.g. Fig. ^ as a function of the incidence 
angle 9, p|| = ppsinO in the region ~ l/{ppa), a be- 
ing typical interlayer distance. To come to the trajectory- 
like picture, one performs averaging in a region of angles 
~ A9. It is the coarse-grain current which should be 
compared with the current density found from the qua- 
siclassical technique. 



II. QUASICLASSICAL CURRENT 

In |t|hc quasiclassical technique,, the current density j 
readso: 



(2.1) 



where e is energy and n a unit vector which shows the 
direction of the momentum, n{e) is the distribution func- 
tion, and the angular and energy resolved current jn(£) 
is found as 



j„(£) =WFnRe(g^(£))^^ 



where vf is the Fermi velocity and is tke retarded 
quasiclassical Green's function (see e.g. Ref.Eil for nota- 
tion). For a given energy e and a parallel momentum 
the X component of the current along the x axis is 
sum of the contributions jni + Jui, , where (ni)^; = cosG 
and (ni/)a; = — cos0 with G being the angle of the 
trajectory defined as sin0 = p \\ / pF- In order to com- 



pare the curre nt w ith Eq. ( 1.28| ), we change the inte- 
gration in Eq. ( ^.l[) as follows: / dfln ^ J dp\\/pF and 
de 2 Jq°° de. The corresponding partial current 
then reads 

PF 



(2.4) 



oo oo 

de [ (e)(l-2n(£)) 

J PF " 




where the Green's functions are taken at the same space 
point, e.g. at the interface shown on Fig.^. 

To evaluate the Green's function g (= g^), we use 
the method of Ref. and express the 2x2 matrix 

g via two-component "wave functions" on classical 
trajectories: 



- 1 



(2.5) 



where = —i(j>'^Ty. These amplitudes-abey Andreev-like 
equation on classical trajectories (seetSEJ)). The index 
± denotes solutions with different asymptotic behaviour: 
the amplitude 0+ — > for the trajectory coordinate x 
going to -|-infinity: a; — > cxd, and 0_ — *■ for x ^oo. 

The Andreev equation needs a boundary condition 
when the trajectory hits an interface and ballistic pieces 
of trajectories are tied by a "knot", see Fig.||a. The 
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FIG. 4: (a) (b) 

The multi-layer structure consisting of A*" superconducting 
layers surrounded by 2 half-infinite superconductors. The or- 
der parameters are Ao,- • • , Ajv+i, the interface transparen- 
cies are Ti, • • • , Tjv+i, and the thicknesses of the layers 
ai, • ■ • , ajv. (a) Scattering on the i'th interface with the trans- 
parency Ti. The angle of the trajectory is O. (b) The zig-zag 
trajectory inside of i'th superconducting layer with the order 
parameter A; . Due to the translational symmetry are certain 
parts of the trajectory equivalent, e.g. 1 ^ 1, 2 ^ 2, • • •. 



(2.2) boundary condition can be formulated via the transfer 



matrix 



M, 



(2.6) 



which relates the wave function on the incoming trajec- 
tory j to that on the outgoing trajectory j'; here and 
below, J denote incoming channels, whereas j' is the out- 
going reflected one, e.g. j = li and j' — 1[ on Fig^a. 
The "knot" transfer matrix Mj'<_j is expressed in terms 
of "across knot" Green's function g^'^k which depends on 
the amplitudes 4)± in the channels on the other side of 
the interfacellil 



1 



T 



2r* 
2 

^'?'(2^) + 0(2.) 
<^(2i)-'?^(2i)+ ' 



1 + R 

- 1 



.92' .2; 



(2.7) 

(2.8) 
(2.9) 



are the transmission and the re- 



52^.2. = 

N = 

where T and R — 
flection probabilities. The indices li and 2i refer to the 
channels on the other side of an interface, see Fig.^a. 
According to Eq.(2.8) the "plus" amplitude 0+ is needed 
only in the outgoing channels and the "minus" amplitude 
in the incoming ones. 
In order to find the Green's function in one of the ex- 
ternal channels (channels which lead to the infinity) we 
need to calculate the amplitudes ^tnd •Pij')- i^i the 

channels inside of the structure. To find these amplitudes 
it is convenient to introduce the total transfer matrix M . 
If one considers the periodic zig-zag trajectory inside of 
i'th layer, see Fig.^ the M is defined as the operator 
connecting the corresponding parts, for example the tra- 
jectories (2i) and (2j) or (2'^) and (2'^). The total transfer 
matrices in the i'th layer shown on Fig.^ read 



(2.10) 



6 



50 



1 COS ■ 



vp COS 9 

e -A, 
A* -e 



Wo sm 



Vp cosQ 



(2.11) 
(2.12) 

(2.13) 



where is the propagator across the i'th layeiS, ^ 



|A,|^Im^ > and 



is the "across knot" 



transfer matrix. The transfer matrices A^ii+i^ii^i and 



Adfi are found in the same way. As shown in 

Ref.0 the quasiclassical Green's function is found as 

9j = - ■ (2-14) 



For example if we want to calculate Green's function in 
the channel (2^) we first find the total transfer matrix 
A^2i^2i connecting ^'s in channels (2^) and (2^) and then 
we calculate g(^2i) from Eg. ( 2. 14 ). 

When the quasiclassical Green's functions 
5(2,)' 5(2'.), 5(1, +i) and g(v^^) are known in each layer, 
one can invert Eq.( |2.5| ) and calculate the amplitudes 
(U+i)- and 



■(20- 



(2;)+' 



1 - fe')22 
(5j')21 

-(5j)l2 
1 - (.9j )22 



(2.15) 
(2.16) 



Using Eqs.(p^^l6|) we can write down the following it- 
erative procedure: 

1. Set all the knot values of the amplitudes 0- in the 
incoming channels and in the outgoing channels 
to the bulk values 



(i')+ 



(1.)- 



(20- 




(2.17) 
(2.18) 



where Ai is the order parameter in the i'th layer, 
andef = e'- |A,|2. 

2. Using the values of the amplitudes (j)±, the "across- 
knot" Green's functions g2'»2i and gi' mii^i are 



constructed from Eq.(2 



3. The "across-knot" Green's functions are substi- 
tuted into Eq.(p.7|) to calculate the knot transfer 
matrices Mji^j, where here and below index j stays 
for 1,; or 2i in i'th layer. 

4. The total transfer matrices -Mj^j and Aip^jt are 



calculated using Eqs.(2.10-2.11) in all layers 
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FIG. 5: Angle resolved current in SISIS structure averaged 
over small range of d6 = 0.09 as a function of energy. The 
angle of incidence 6 = 7r/4. The solid and dashed lines show 
for the quasiclassical and coarse-grain BdG current, respec- 
tively. The bars show the fluctuation of the high-resolution 
BdG current around its average value. The order parameter 
in the leftmost superconductor is Aq. The order parameters 
in the next layers read Ai = e"^''*Ao, A2 = iAo. The inter- 
face transparencies Ti = 0.1, r2 = 0.5. the thickness of the 
layer a = vf/^o- 



5. Using the values of the total transfer matrices from 
step ^. the quasiclassic al Gr een's functions gj and 
gj' are found from Eq.( ^.14 ). 



6. At this stage the quasiclassical Green's function gj 
and gji are known in all layers and the amplitudes 
^ (i)^ and (1)^')+ can be evaluated using Eq.(2.15- 



2.16) 



7. Continue from point ^ until the convergence is 
reached. 

After the last iteration the internal amplitudes 4'(j)- 



and 



(i')+ 



are known and one constructs the "knot" 



transfer matrix Mi/^i on the leftmost external inter- 



face. Since </>(!')+ and (t>n \- are the bulk superconductor 



amplitudes (see Eqs.( 2.17 - 2.18| )) one calculates also the 
Gre en's functions 5(i), 5(i') and the partial current in 
Eq.(|l. 



III. RESULTS 

In this section we present results of the calculations 
for typical parameters of the multi-layer structure such 
as distribution of the order parameter, thicknesses and 
number layers (barriers), and the strength of the barriers. 
Thicknesses of the layers are usually of order of vf/Aq 
where Aq is the gap in the leftmost layer. The value 
Pf — W^Aq/vf is chosen for the Fermi momentum. 

A typical high resolution angular dependence of the 
current jBdcid) has been already shown in Fig. ||. As 
expected, it does not have any resemblance to the smooth 
quasiclassical behaviour. However, after the coarse grain 
averaging, i.e. on a low resolution level, jBdG is in perfect 
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FIG. 6: Angle resolved current in a three-barrier structure 
SISISIS as a function of angle 6 as found from the BdG 
equation. The order parameter in the leftmost supercon- 
ductor is Aq. The order parameters in the next layers read 
Ai = iAo,A2 = — Ao,A3 = — iAq. The thicknesses of the 
two internal layers a\ = UF/Ao,a2 = vf/Aq. The energy 
£ — I.IAq. The interface transparencies Ti — 0.1, r2 — 
0.5, Ta — 0.5. The solid line corresponds to the quasiclas- 
sical current, a constant in this narrow interval. The dashed 
line shows the BdG current. 



agreement with quasiclassics, see Fig.| The "error-bars" 
in Fig.^ show the mean square fluctuation of the high 
resolution current around its coarse-grain average. In the 
double-barrier case, the agreement exists for any angle 
9 and energy e, and for any set of parameters of the 
structure, A's and the barrier's strength. 

Contrary, noticeable deviations from the quasiclassical 
solutions are seen when there are more than two barriers. 
Here we present results only for three-barrier structures, 
more complicated systems show qualitatively same fea- 
tures. 

For a three-interface structure, a typical high resolu- 
tion angle dependence of JbcIG is shown in Fig. ^ In 
Fig.^ we plot coarse-grain BdG-current together with 
the quasiclassical curve for slightly different geometries: 
Fig.^(a) refers to a symmetric case when the two inter- 
nal layers have exactly the same thickness whereas in 
0(b) the thicknesses are 10 percent different from each 
other. In the both cases, one sees a clear deviation of the 
quasiclassical curve from the "exact" one, the deviation 
lesser in asymmetric geometry. 

For the same geometries and the order parameters, the 
energy dependence of the current integrated with respect 
to the incident angle are shown in FigJ^. Disagreement 
between the exact and quasiclassical results are clearly 
seen, again stronger in the symmetric case. 



IV. DISCUSSION 

The results of the previous section clearly show that in 
some geometries the quasiclassical theory does not repro- 
duce the "exact" results derived from the BdG equation. 
The two approaches agree only qualitatively. The quan- 





FIG. 7: (a) (b) 

Angle resolved current averaged over small range of = 0.09 
as a function of energy. The angle of incidence 6 — 7r/4. 
The solid line stands for the quasiclassics and the dashed line 
for the BdG current. The bars show how the BdG current 
changes around its average value. The order parameter in 
the leftmost superconductor is Aq. The order parameters in 
the next layers read Ai = iAq, A2 = — Aq, A3 = — jAq. The 
interface transparencies Ti = 0.1, T2 = 0.5, T3 = 0.1. (a) 
Symmetric case: the thicknesses of the two internal layers 
ai — «F/Ao,a2 = vf/Ao- (b) Non-symmetric case: ai = 
0.9 VF/Ao,a2 = 1.1 vf/Ao. 
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FIG. 8: (a) (b) 

The total current as a function of energy. The solid line stands 
for the quasiclassical and the dashed line for the BdG cur- 
rent. The order parameter in the leftmost superconductor 
is Aq. The order parameters in the next layers read Ai = 
iAo,A2 — — Ao,A3 — —iAo- The interface transparencies 
Ti = 0.1, T2 = 0.5, T3 = 0.1. (a) Symmetric case: the thick- 
nesses of the two internal layers ai = VF/Ao,a2 — vf/Aq. 
(b) Asymmetric case: ai = 0.9 VF/Ao,a2 — 1.1 vf/Aq. 



titative discrepancy much exceeds the corrections to the 
quasiclassical theory of order of l/ppa ~ ^/pfVf which 
one might expect. Below, we present our understanding 
of physics behind the discre p anc y . 

As in our earlier papersliJIia'tZl, we ascribe the failure 
of the quasiclassical theory to the presence of interfering 
paths or, in other words, loop-like trajectories. From this 
point, the validity of the quasiclassical theory in the two- 
barrier case (see Fig]|) is in accordance with our expecta- 
tions. Indeed, in this simple geometry, the classical path 
shown in Fig. |^a, is effectively one diruansional (tree- 
like trajectory in the terminology of Ref.liJ) in the sense 
that there is only op£ path connecting any two points. 
As discussed in Ref.Eil, one is then able to factorize the 
full propagator G{x,x'), x x' labeUing the points on the 
tree-like trajectory, as G{x,x') = g{x,x')&cp[ipF^xx'], 
where Cxx' is the length of the path along the tree-like 
trajectory connecting x and x' , and g{x,x') is a slowly 
varying quasiclassical (2-point) Green's function. In this 
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to the full Andreev reflection amplitude Ae is evaluated 
in appendix 



FIG. 9: (a) (b) 

Classical trajectories. The trajectory is built of ballistic pieces 
"tied" by scattering on the two interfaces (knots). The arrows 
show the direction of the momentum. There are no interfering 
paths in a double-layer case (a). Loops exist in a three-barrier 
system (b). 
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FIG. 10: (a) (b) 

Simplest loop like paths contributing to the Andreev reflec- 
tion amplitude in a symmetric double-barrier structure. The 
incoming electron comes from the left and is reflected as a 
hole after going around the loop anti-clockwise (a) and clock- 
wise (b). Arrows on the continuous lines show the direction 
of momentum. Letters "e" and "h" show the character of the 
excitation - electron or hole. The arrows under the letters 
show the direction of propagation. 



case, derivation of the quasiclassical equation does not 
meet any difSculty, and the theory gives valid results. 

Turning now to the symmetric three-barrier structure, 
the trajectory is not simple tree- like (see Fig.^(b)): there 
are loops and, therefore, interfering paths. Then, one 
may expect corrections to the quasiclassical theorVj-which 
are not controlled by the quasiclassical parametertZI. The 
origin of the corrections and its relation to the existence 
of loops, can be understood from the following semi- 
quantitative arguments. 



In the BdG-picture the current Eq. (1.28) is expressed 
via the element of the S-matrix that is the Andreev 
reflection amplitude = Syi- Scattering amplitudes 
can be presented as a sum of partial amplitudes, each of 
which corresponds to a particular path of the particle. 
Among others, there are closed paths shown in Fig.[ic|. 
The contribution A^°°p*' of the paths in FigJlO|(a) and (b) 



^loops ^J^^ 



2i(cos Qpp{^a\—a'2)^ 



(4.1) 



where is a coefflcient defined in Eq. (B5), a\ (0,2) is 
the distance from the barrier 1 to barrier 2 (from 2 to 
3); r\ and ra are the amplitudes of reflection from the 
barrier 1 and 3 respectively. Note that this simplest loop 
survives the coarse grain averaging only if a\ « a-i. We 
understand the larger deviation from quasiclassics seen in 
the symmetric case compared with an asymmetric one, 
as due to the contribution of the simple loop. 

The existence and importance of this contribution can 
be checked exploiting the fact that it is sensitive to the 
phase of the reflection coefficients and the length of the 
path on the scale of 1/pf- In Fig.[lT|, we show the cur- 
rent for different signs of the barrier strength A3 = ijAaj, 
changing the phase of = —i\^/{pFx + iAsV) but leav- 
ing the reflection probability intact. In Fig.O, we plot 
the change of the exact and quasiclassical currents upon 
tiny variation of the right layer thickness (correspond- 
ing to TT-change of the Fermi phase factor in Eq. ( [4.l|) ). 
While quasiclassical current remains intact, clearly seen 
changes are observed in the exact curr ent with the order 
of magnitude consistent with Eq. (4.1). 

To avoid confusion, we remind that we deal with 
coarse-grain averaged currents, and therefore the ob- 
served sensitivity to the thickness and the phase of reflec- 
tion has nothing to do with the size effects (due to the 
commensurability of the thickness and the Fermi wave 
length) well known in the normal case. We note also, 
that the loops in Fig.|l^ do not exist in the normal state 
because the electron-hole conversion on the interface 2 
would not be possible. 



We assert that the loop contribution Eq. ( [4.1| ) is chiefly 
responsible for the deviations from the quasiclassical the- 
ory. Obviously, this contribution cannot be grasped by 
quasiclassics since A^°°^^ is sensitive to the phase of the 
reflection amplitudes ri anf ^^a, whereas the quasiclassi- 
cal boundary conditionO'E3'L3 contain only the probabil- 
ities |rp and This is our argument supporting our 
interpretation of the numerical results. Note that the in- 
terpretation is consistence with the observation that the 
deviation from quasiclassics are significantly smaller in 
the asymmetric case Fig||: There, the simple loops are 
absent and the deviations from the conventional quasi- 
classics come from higher order loops (like that analysed 

iiM). 



V. CONCLUSIONS 

In this paper we have examined the applicability of 
quasiclassical theory for description of multiple interface 
scattering by comparing quasiclassical solutions with "ex- 
act" ones, extracted from coarse-grain averaged solutions 
to the Bogoliubov - de Gennes equation. We see that the 
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Quasiclassics 




FIG. 11: Angle resolved current averaged over small range 
oi d9 — 0.09 as a function of energy. The angle of inci- 
dence 9 = 7r/4. The solid line stands for the quasiclassi- 
cal current. The other two lines correspond to A3 positive 
or negative. The order parameter in the leftmost supercon- 
ductor is Ao. The order parameters in the next layers read 
Ai — iAo, A2 — — Ao, A3 = — iAq. The interface transparen- 
cies ri=0.5,r2 = 0.5,r3 = 0.5. The thicknesses of the two 
internal layers ai — vf/Ao,o-2 = vp/^o- 
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FIG. 12: The differences between currents in SSSS setups dif- 
fering only by the thickness of the last layer. The thicknesses 
are ai = a2 — wf/Aq in the first case and ai = vf/Aq and 
a2 = 1.002wf/Ao in the second case. The order parameter in 
the leftmost superconductor is Aq. The order parameters in 
the next layers read Ai = lAo, A2 = — Ao, A3 = —iAo. The 
interface transparencies Ti = 0.5, r2 = 0.5, r3 — 0.5. 



two approaches agree in simple geometries (one or two 
interfaces) but show noticeable discrepancy when three 
or more interfaces are present. This gives an example of 
a physical system where quasiclassical technique fails to 
give quantitative description with its expected accuracy 

1/pf£,o- As we understand it, the failure of the qua- 
siclassical theory occur when classical trajectories form 
closed loops (interfering paths) after sequential reflection 
and transmission accompanied by electron-hole conver- 
sion. Th|is gives additional support to the point of view 
of Ref.EJ that the derivation of the quasiclassical tech- 
nique is possible only under the assumption of a simply 
connected topology, tree-like, of classical trajectories. 

The main goal of this paper has been to demonstrate 



the existence of noticeable deviations from the quasiclas- 
sical theory in conditions where one might expect it to 
give fully reliable results. For this purpose we have cho- 
sen the simplest "exact" method, the Bogoliubov - de 
Gennes approach, where the superconductivity enters via 
the mean-field order parameter A(r) and scattering due 
to either impurity or surface roughness is not included. 
It is now time to discuss to what extent our results are 
sensitive to the simplifications. 

The truly "exact" theory of (phonon-mediated) super- 
conductivity, for which the quasiclassical technique is an 
approximation, is the set of Gor'kov-Eliashjberg equations 
for disorder averaged Green's function QJE^. The equa- 
tion of motion, (e — Tig) Qe — 1; contains the operator 
He which has the same structure as the Bogoliubov - de 
Gennes Hamiltonian Eq. (1.1), with the order parameter 
replaced by the anomalous part of the electron-phonon 
self-energy. Apart from the e-dependence of He , the only 
qualitative difference is that the Ti.^ has a non-Hermitian 
part coming from the self-energies. The latter accounts 
for the impurity and electron-phonon scattering and, cor- 
respondingly, for a finite life time of the excitation with 
a given momentum, T(e) - the scattering-out time. Due 
to the similarity of the above operators, the interference 
contribution to observables obtained with the help of the 
Green's function technique or simple minded mean-field 
Bogoliubov - de Gennes equation, would give same re- 
sults. A finite life time is the only important feature 
missing in the Bogoliubov - de Gennes approach, and 
below we discuss its role. 

Clearly, the interference of waves having travelled dif- 
ferent paths occurs if only the decay length is not too 
small compared with the path lengths. By virtue of the 
optical theorem, the waves corresponding to the ballistic 
trajectories decays on the distance ~ vpT. In practice, 
T is controlled by the bulk impurity scattering, so that 
the loops in Fig. |l^ may contribute only if the interlayer 
distances are less or of order of the impurity mean path. 
As far as interface imperfection is concerned, short-range 
surface roughness is expected to play the role similar to 
that of the bulk impurity scattering: Although the inter- 
face reflection becomes partially diffusive, the coherent 
specular component, with the irttepaity proportional to 
the Fuchs's parameter V, is finitec^lHJ. Thus, the picture 
of trajectories as that sho wn in Fig. ^ remains mean- 
ingful as well as Eq. ( |4.lD if corresponding attenuation 
factors are inserted. We see that although microscopic 
roughness and bulk scattering suppress the interference 
effect, it survives disorder averaging if the disorder is not 
too strong. 

When a long-range roughness is present, that is the 
layer thicknesses are slowly vary ing, the global value of 
the loop contribution in Eq. (4.1) averages to zero (if the 
thickness modulation exceeds Xp)- However, the inter- 
ference of paths forming the loops is sensitive mainly to 
the local value of the thicknesses, so that the interference 
effects can be seen in the spatial fluctuations of the local 
current, which is a measurable quantity. In the quasiclas- 
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sical theory, the roughness would reveal itself only of the 
thickness modulation is somehow comparable with the 
coherence length. Since the loops are sensitive to varia- 
tions of geometry on the scale of Xp, their contribution 
to the fluctuations is expected to be much larger than 
that in the quasiclassical theory. 

We see that the deviation from the quasiclassical the- 
ory due to the interference effects, although most pro- 
nounced in the idealized model exploited in the paper, are 
observable in realistic conditions if the disorder is not too 
strong. After disorder averaging, the interference contri- 
bution is observable if the mean free path exceeds the 
interlayer distance and the rough interfaces have not too 
small coefficient of the specular reflection. Detail analy- 
sis of the loop contribution to mesoscopic fluctuations is 
beyond the scope of the paper. 

Apart from disorder, the interference contribution may 
be suppressed by energy integration. Indeed, the inte- 
gration with the Fermi-Dirac distribution function corre- 
sponding to the temperature T is equivalent to the Mat- 
subara summation, that is energy variable e assumes dis- 
crete imaginary values, multiples of inT. Then, the waves 



decay on the length ^ ^q, and, consequently, the loops 
contribute to equilibrium properties only if the paths are 
shorter than the coherence length. 

Summarizing, we have shown that the quasiclassical 
technique fails in geometries where classical trajectories 
form closed loops. In particular, this happens when it is 
applied to the Josephson multi-layer structure with num- 
ber of semi-transparent interfaces larger than 2. This 
conclusion does not undermine the conventional quasi- 
classical technique but only limits its applicability in 
some special geometries where the interference of classi- 
cal paths cannot be neglected. Within the quasiclassical 
approach, the interference can be incorporated into the 
theory with the help of the method suggested in Ref. |l^ . 
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APPENDIX A: THE TRANSFER AND 
SCATTERING MATRICES 

The Bogoliubov - de Gennes equation is a second-order 
differential equation. In the one-dimensional case consid- 
ered here, it can be reduced to the first order equation 
for an "extended" wave functions ^' which is built of the 



where a, b, c, and d are real 
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wave function -0 and its derivative ptp, 



( ^ 



V 



dx 



(Al) 



Since ■0 has two components, the extended wave function 
is a 4-coniponent column. 

In ter ms of 5*, the quasiparticle current density 
Eq. (1.6) reads 



3 



qp 



1 

2m 



(A2) 



where is the PauH matrix operating in the space ip 
pip, and, as before, acts in t he u — v space. 
The charge current Eq. (E3) can be found as 



e 2m 



(A3) 



The extend ed wa ve function corresponding to the plane 
waves in Eq. (1.11) have the form 



*i.<T(a;) 



„icrp^x ^ 



(A4) 



where the 4-component amplitudes $^cr may be taken as 



1 



or in a more concise form, 
where 



1 



/2p. 



1 

<7Pu 



(A5) 



(A6) 



(A7) 



The amplitudes '^^a are normalized to the probability 
flux Eq. (A2) equal to l/2m (for E > |A|); accordingly, 
the plane waves in Eq. (A4) are normalised to the S- 
function of energy. 

The conjugated amplitudes defined as 



*L = V'i «) 0* „ where 



'^'P- V 2 V'p. 



(A8) 

{(j)]y p = a(j)^ pdx) satisfy the following orthogonality and 
completeness relations 



(A9) 



Due to the orthogonality property, <i>^ ^ projects a gen- 
eral superposition ^'(a;) to the plane wave ^'^ ^•(a;)- Using 
this argument, one constructs the evolution operator Ua, 
a 4 X 4 matrix, which relates the wave functions at the 
points X and x + a, ^>{x) = Ua'^ix + a). 



(AlO) 



In the model where the potential barrier V{x) separat- 
ing the layers is a ^— function, V{x) = —6{x), the two- 
component wave function ■0(x) is continuous at a; = 0, 
whereas the derivatives suffer the jump: V'loi = Xip{0). 
In terms of the extended wave function 5* Eq. (Al), the 
interface matching condition reads 



V 



1 
2iX 1 



(All) 



It is implied in Eq. (All) that each matrix element of V 
is multiplied by the unit matrix in the u — v space so that 
V is actually a 4 x 4EI. 

The transfer matrix, M, relates the extended wave 
function, on the left side of the multi-layer structure 
to that on the right side 



It is given by the ordered product, 

M ^ ViUi 2X>2 ■ • -Un-i nT>n 



(A12) 



(A13) 



of the matrices Vk Eq. ( All ) corresponding to the fc— th 
interface, k — 1, . . . N, and the evolution matrices Uk,k+i 
accounting for the propagation from the k + 1-th to fc-th 
potential barrier. 

The elements of the S-matrix can be found via the 
transfer matrix. For this, we take advantage of the com- 
pleteness relation in Eq. (|A9[) and present M. as 



M 



(A14) 



where we denote /x the set {v,(j) and /x' — {v' ,a'). The 
elements of the matrix M read 



'mm' 



(A15) 



where again fj. stands for {v, a) and '^^a and ^Vo- are 
the 4-component amplitudes Eq. ( |A5| ) and Eq. ([A^), re- 
spectively. The meaning of the M-matrix is that it is 
the transfer matrix in the plane wave representation: 



(fl) 



where C's are the coefficients in 



the expansion = ^ c\t 



PresentiM 
block formc^l. 



found from Eqs. ( A13| ), and (A15), in a 



M 



A D 
B C 



(A16) 



the S'-matrix expressed via 2x2 matrices A,B,C, and 
D, reads 



S = 



C-BA-^D 
-A-^D 



(A17) 



The matrix element Skn gives the amplitude of the scat- 
tering from the n— th incoming state listed in Eq. (1.1£) 
to the k-ih final state in Eq. (1.20). 
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APPENDIX B: LOOP CONTRIBUTION 

To find contribution of the loops, we first analyze the 
elementary process that is the scattering on an isolated 
barrier. Consider a single barrier in between two semi- 
infinite homogeniuos superconductors left (L), and right 
(R). The barrier is characterized by the reflection r and 
transmission amplitudes t. As required by unitarity, rt* + 
*< = and i? + r = 1 where R = H^T = The free 



wave functions arc listed in Eqs. (1.19), and (1.20). 

The scattering matrix calculated by the method de- 
scribed in in Appendix Q reads 



S = 7 



t(3 
\ r*tal3* 



where 



a 



Itpa* t(3* r*ta*l3 

r*P(3* t*ra(3* t*(3 
t*ra*l3 r/3/3* \t\'^a 

t*l3* a*|t|2 r*Pl3* 



(Bl) 



1 - 



Electron (^-^^"^^ 



(B2) 

and hole (?/;L^''^'') wave functions are 



defined in Eq. (1.14). By their physical meaning, a is 
the amplitude of the Andreev reflection on the SS' in- 
terface in the absence of barrier, and /3 is the transmis- 
sion amplitudes. Outside the energy gap, a(3* + a* (3 — 
, aa*-\-l3f3* = 1 as required by the quasiparticle current 
conservation. As no surprise, the S-matrix in Eq. (Bl ) has 
the same structure as that derived in Ref.|| for the NIS 
interface. 

The loop contribution to the Andreev reflection A^°°^^ 



is the sum processes shown 

in Fig.^O|(a) and (b). The amplitudes of each of the pro- 
cesses is the product of factors accumulated along the 
path. The rules to find the factors are as follows: 



(i) The factor which corresponds to the ballistic part 
of the trajectory is ey:jp[ivp^\x f — Xi \\ wh ere p^, is the 
X— component of the momentum Eq. (1.12), v — ±is the 
type of the excitation (electron, or hole, "— " ), and 
Xi and X/ is the initial and final value of the x-coordinate. 
One can prove this formula taking into consideration that 
the electron propagates in the direction of momentum, 
and the hole in the opposite direction. (The phase ac- 
cumulated due to a displacement in the p| [-direction can 
be omitted since p\ \ is the same for all the ballistic pieces 
and the path is closed.) 

(ii) For an interface scattering event, the factor is the 
element of the S-matrix corresponding to the initial and 



final states in Eq. (Bl 



Looking at Fig. |10|(a,b) and using these rules, one gets 

/l(a) _ c(l) „jp+ai 0(2) _-ip-a2 c(3) -ip_a2 c(2) „ip+ai c(l) 
— '-'23 ^ '-'14 ^ '-'22 ^ '-'41 ^ '-'31 

where superscript in S'^^\ k = 1,2, 3, labels the interface 
(see Fig.pX]|) , and ai (02) is the distance from the barrier 1 
to barrier2 (from 2 to 3); the momentum pj- is calculated 
for the parameters of the corresponding layer. 

Finally, substituting the elements of the S-matrix 
Eq. (Bl), one gets. 



^loops ^ ^(a) ^ ^(b) 



(B3) 



where the coefficient A reads 

^ = -ai7h273|r-2iii2a2/3i/32/33pe'^'^(«^'^^+«^'^^' . 

(B5) 

Here, 9 is the angle between the direction of th e tra jec- 
tory and the x— axis, and ^1,2 is defined by Eq. ( |l.l3| ) for 
the corresponding layer. 



